! This is a 1D advection example using square initial condition and periodic
! boundary condition for BiMocq advection scheme.
!
! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2024-03-21: Initial creation.

program bimocq_adv_1d_case

  use adv_1d_test_case_mod
  use bimocq_common_mod

  implicit none

  real, allocatable :: xmap(:)            ! Characteristic mapping
  real, allocatable :: xflx(:)            ! Flux of characteristic mapping
  real, allocatable :: rho0(:)            ! Initial tracer density
  real, allocatable :: rho (:)            ! Tracer density being advected at cell centers
  real x0, w1, w2
  integer, parameter :: ns = 3            ! Stencil width
  integer i, i0
  character(30), parameter :: scheme = 'bimocq_1d'

  namelist /params/ nx, nt, dt, u

  call get_command_argument(1, namelist_path)
  inquire(file=namelist_path, exist=is_exist)
  if (is_exist) then
    open(10, file=namelist_path)
    read(10, nml=params)
    close(10)
  end if

  allocate(xmap(1-ns:nx+ns))
  allocate(xflx(1-ns:nx+ns))
  allocate(rho0(1-ns:nx+ns))
  allocate(rho (1-ns:nx+ns))

  call adv_1d_test_case_init('square', ns, rho)
  call reset()
  call output(scheme, 0, ns, nx, x, rho)

  ! Run integration.
  print *, time_step, sum(rho(1:nx))
  do while (time_step < nt)
    call advect_map()
    call interp_rho()
    call advance_time()
    if (mod(time_step, 20) == 0) call reset()
    call output(scheme, time_step, ns, nx, x, rho)
    print *, time_step, sum(rho(1:nx))
  end do

  call adv_1d_test_case_final()

  deallocate(xmap, xflx, rho0, rho)

contains

  subroutine advect_map()

    call upwind(xmap, xflx)
    do i = 1, nx
      xmap(i) = xmap(i) - dt / dx * (xflx(i) - xflx(i-1))
    end do
    call apply_bc(ns, nx, xmap)
    xmap(1-ns:0    ) = xmap(1-ns:0    ) - xrange
    xmap(nx+1:nx+ns) = xmap(nx+1:nx+ns) + xrange

  end subroutine advect_map

  subroutine interp_rho()

    do i = 1, nx
      x0 = limit_coord(xmin, xmax, xrange, xmap(i))
      call interp1d(1, ns, nx, x, rho0, x0, rho(i))
    end do

  end subroutine interp_rho

  subroutine reset()

    do i = 1, nx
      xmap(i) = x(i)
    end do
    call apply_bc(ns, nx, xmap)
    xmap(1-ns:0    ) = xmap(1-ns:0    ) - xrange
    xmap(nx+1:nx+ns) = xmap(nx+1:nx+ns) + xrange
    if (time_step == 0) then
      rho0 = rho
    else
      rho0(1:nx) = rho(1:nx) / sum(rho(1:nx)) * sum(rho0(1:nx))
    end if
    call apply_bc(ns, nx, rho0)

  end subroutine reset

  subroutine upwind(q, f)

    real, intent(in ) :: q(1-ns:nx+ns)
    real, intent(out) :: f(1-ns:nx+ns)

    real, parameter :: beta = 0.5

    do i = 0, nx
      f(i) = 0.5 * (u * (q(i+1) + q(i)) - beta * abs(u) * (q(i+1) - q(i)))
    end do

  end subroutine upwind

end program bimocq_adv_1d_case
